Dla wszystkich wymienionych powyżej sformułowań wariacyjnych w celu rozwiązania ich metodą elemntów skończonych konieczne jest wygenerowanie układu równań liniowych. W szczególności konieczne jest wygenerowanie całek do układu równań liniowych. Pamiętając, że wiersze i kolumny w globalnym układzie równań odpowiadają poszczególnym funkcjom bazowym, w naszym przykładzie odpowiadają one funkcjom B-spline, generujemy następujący układ równań
\( \begin{bmatrix} & \color{red}{B^x_{1,p} B^y_{1,p}} & \cdots &\color{red}{ B^x_{N_x,p}B^y_{1,p}} \\ \color{blue}{ B^x_{1,p}B^y_{1,p}}& a(\color{blue}{ B^x_{1,p}B^y_{1,p}}, \color{red}{ B^x_{1,p}B^y_{1,p}}) & \cdots & a(\color{blue}{ B^x_{1,p}B^y_{1,p}},\color{red}{ B^x_{N_x,p}B^y_{N_y,p}}) \\ \color{blue}{ B^x_{2,p}B^y_{1,p}} & a(\color{blue}{ B^x_{2,p}B^y_{1,p}},\color{red}{ B^x_{1,p}B^y_{1,p}}) & \cdots & a(\color{blue}{ B^x_{2,p}B^y_{1,p}},\color{red}{ B^x_{N_x,p}B^y_{N_y,p}}) \\ \vdots & \vdots & \vdots & \vdots \\ \color{blue}{ B^x_{N_x,p}B^y_{N_y,p} } & a(\color{blue}{ B^x_{N_x,p}B^y_{N_y,p}}, \color{red}{ B^x_{1,p}B^y_{1,p}}) & \cdots & a(\color{blue}{ B^x_{N_x,p}B^y_{N_y,p}},\color{red}{ B^x_{N_x,p}B^y_{N_y,p}}) \\ \end{bmatrix} \begin{bmatrix} u_{1,1} \\ u_{2,1} \\ \vdots \\ u_{N_x,N_y } \\ \end{bmatrix} \)
\( = \begin{bmatrix} l(B^x_1(x)*B^y_1(y)) \\ l(B^x_1(x)*B^y_2(y)) \\ \vdots \\ l(B^x_{N_x}(x)*B^y_{N_y}(y)) \\ \end{bmatrix} \)
Funkcje B-spline które nie mają wspólnej dziedziny (które znajdują się daleko od siebie na siatce) generują zerowe całki, czyli zerowa wartości macierzy. Podczas generacji macierzy nie należy więc iterować po wierszach i kolumnach macierzy, za to należy iterować po elementach siatki obliczeniowej, oraz po funkcjach bazowych które są określone na danym elemencie. Wynikowe fragmenty całek należy dopisywać do globalnej macierzy. Wówczas to nasza generacja całek omijać będzie wszystkie niezerowe wyrazy macierzy.
Algorytm zawiera sześć zagnieżdżonych pętli:
1 for \( nex=1,N_x \) (pętla po elementach wzdłuż osi \( x \))
2 for \( ney=1,N_y \) (pętla po elementach wzdłuż osi \( y \))
3 for \( ibx1=1,p+1 \) (pętla po \( p+1 \) B-spline elementu wzdłuż osi \( x \))
4 for \( iby1=1,p+1 \) ( pętla po \( p+1 \) B-spline elementu wzdłuż osi \( y \))
5 \( i = f(nex,ibx1) \) ( oblicz indeks B-spline wzdłuż osi \( x \))
6 \( j = f(ney,iby1) \) ( oblicz indeks B-spline wzdłuż osi \( y \))
7 \( irow = g(nex,ibx1,ney,iby1) \)(globalny indeks 2D B-spline'a)
8 \( l(irow)+= l(B^x_{i,p}B^y_{j,p}) \)
9 for \( ibx2=1,p+1 \) ( pętla po \( p+1 \) B-spline elementu wzdłuż osi \( x \))
10 for \( iby2=1,p+1 \) (pętla po \( p+1 \) B-spline elementu wzdłuż osi \( y \))
11 \( k = f(nex,ibx2) \) (oblicz indeks B-spline wzdłuż osi \( x \))
12 \( l = f(ney,iby2) \) (oblicz indeks B-spline wzdłuż osi \( y \))
13 \( irow = g(nex,ibx2,ney,iby2) \) (globalny indeks 2D B-spline'a)
14 \( M(irow,icol)+= a(B^x_{i,p}B^y_{j,p},B^x_{k,p}B^y_{l,p}) \)
Funkcje \( f(nex,ibx) \) obliczają indeks
funkcji \( ibx \)-tej B-spline rozpiętej na elemencie wzdłuż osi \( x \) numer \( nex \). Dla jednorodnych B-spline'ów stopnia \( p \) indeks ten wynosi \( nex+ibx-1 \).
Funkcje \( g(nex,ibx,ney,iby) \) obliczają globalny indeks dwuwymiarowego B-spline'a. W przypadku jednorodnych funkcji B-spline wzór ten wynosi \( irow = i*j \) gdzie \( i=f(nex,ibx) \) oraz \( j=f(ney,iby) \).
Z koleji obliczanie całek realizowane jest w ogólnym przypadku przez całkowanie numeryczne z wykorzystaniem punktów i wag kwadratury, na przykład kwadratury Gaussa. Całkę jednowymiarowaz wielomianu na przedziale da się obliczyć dokładnie jako suma wartości tego wielomianu w wybranych punktach kwadrtury rozpostartych na przedziale, przemnożonych przez odpowiednie wagi. Wyprowadzenie i teoria kwadratur numerycznych wykracza poza zakres tego podręcznika.
Nazwa kwadratur Gaussa pochodzi od nazwiska niemieckiego matematyka który wynalazł tą metodę całkowania w 1826 roku.
[1]
Wybór kwadratury zależy od rodzaju funkcji bazowych, i możliwa jest redukcja liczby punktów kwadratury stosując specjalne kwadratury dedykowane poszczególnym funkcjom bazowym. W ogólnym przypadku na pojedynczym elemencie funkcje bazowe są wielomianami i stosowanie kwadratur Gaussa jest uzasadnione.
Poniższy przykład ilustruje tworzenia punktów i wag kwadratury Gaussa na przedziale \( [0,1] \) gwarantującej dokładne policzenie całki z funkcji dwuwymiarowych B-spline stopnia nie większego niż 10.
m_points[1]= (1.0-0.973906528517171720077964)*0.5;
m_points[2]= (1.0-0.8650633666889845107320967)*0.5;
m_points[3]= (1.0-0.6794095682990244062343274)*0.5;
m_points[4]= (1.0-0.4333953941292471907992659)*0.5;
m_points[5]= (1.0-0.1488743389816312108848260)*0.5;
m_points[6]= (1.0+0.1488743389816312108848260)*0.5;
m_points[7]= (1.0+0.4333953941292471907992659)*0.5;
m_points[8]= (1.0+0.6794095682990244062343274)*0.5;
m_points[9]= (1.0+0.8650633666889845107320967)*0.5;
m_points[10]= (1.0+0.9739065285171717200779640)*0.5;
m_weights[1]=0.0666713443086881375935688*0.5;
m_weights[2]=0.1494513491505805931457763*0.5;
m_weights[3]=0.2190863625159820439955349*0.5;
m_weights[4]=0.2692667193099963550912269*0.5;
m_weights[5]=0.2955242247147528701738930*0.5;
m_weights[6]=0.2955242247147528701738930*0.5;
m_weights[7]=0.2692667193099963550912269*0.5;
m_weights[8]=0.2190863625159820439955349*0.5;
m_weights[9]=0.1494513491505805931457763*0.5;
m_weights[10]=0.0666713443086881375935688*0.5;
Kwadratury oryginalnie tablicowane są dla przedziału \( [-1,1] \) tak więc tutaj przeskalowujemy je do przedziału \( [0,1] \) dodając 1 do współrzędnej punktów oraz dzieląc współrzędne punktow przez 2 (czyli mnożąc przez 0.5).
Wówczas całkowanie numeryczne wykonywane w linii 14 (i kolejnych) w powżzszym kodzie wymaga dwóch dodatkowych pętli po punktach i wagach kwadratury Gassua i wygląda ono następująco:
15 for \( m=1,p/2+1 \) (pętla po punktach kwadratury Gaussa wzdłuż osi \( x \) dla B-spline'ów stopnia \( p \))
16 for \( n=1,p/2+1 \) (pętla po punktach kwadratury Gaussa wzdłuż osi \( y \) dla B-spline'ów stopnia \( p \))
17 \( M(irow,icol)+=\\ a(B^x_{i,p}(m\_points[m])B^y_{j,p}(m\_points[n]), B^x_{k,p}(m\_points[m])B^y_{l,p}(m\_points[n])) \\ *m\_weights[m]*m\_weights[n]*area \)
W szczególności dla naszego przykładowego problemu transportu ciepła na obszarze w ksztalcie litery L mamy
\( M(u,v) = \left(\frac{\partial B^x_{i,p} }{\partial x }(m\_points[m])\frac{\partial B^x_{k,p} }{\partial x }(m\_points[m])\right)*m\_weights[m] \\ + \left(\frac{\partial B^y_{j,p} }{\partial y }(m\_points[n]) \frac{\partial B^y_{l,p} }{\partial y }(m\_points[n]) \right)*m\_weights[n] \)
W powyższym wzorze obliczamy wartości pochodnych z jednowymiarowych funkcji B-spline w punktach kwadratury Gaussa \( (m\_points[m],m\_points[n]) \) przemnożone przez wagi kwadratury Gaussa \( *m\_weights[m]*m\_weights[n] \).